FEATURE ARTICLE 



by Michael Dvorsky 



Rooting Around 

Integer Square Roots on Small Processors 

There are numerous ways to compute square roots. Have you ever tried the sum-of-odds 
method? What about the bisection method? Michael goes over the basics and shows you 
how to choose algorithms to suit your particular applications. 



Square roots and. I go way back. I 
was in high school when so-called , 
"four-banger" calculators became 
affordable. These calculators had only 
four standard functions: addition, sub- 
traction, multiplication, and division. 
But what if I wanted to do more exotic 
computations like square root? My - 
local bookstore had several books 
about calculator tricks, so I picked up a 
couple to see what I could do with my 
new toy. (Remember the excitement of 
punching in "07734" and turning the 
calculator upside-down to read the 
"hELLO" message?) One of these books 
described the Newton-Raphson method 
for computing square roots, which is an 
iterative procedure that used only the 
standard four functions. I memorized 
the algorithm, and could still punch it 
out long after I had forgotten the 
tedious high school method for com- 
puting square roots by hand. 

My interest in the process of comput- 
ing square roots was renewed a decade 
ago when I read an exchange on the 
Circuit Cellar BBS. Somebody asked 
about the sum<-of-odds method, which 
is a simple and seemingly fast way to 
compute integer square roots. Circuit 
Cellar columnist Ed Nisley suggested 
the use of bisection (or binary search) to 
quickly find the square root using a 
processor with hardware multiply. 

Over the years, I returned to the 
square root problem on numerous 
occasions. I came up with several algo- 
rithms that I later learned were well 
known in the field. In this article, I'll 
describe how these algorithms work. 



I'll also help you determine which 
algorithm is best for your application. 

WHY ROOT? 

Before getting into the details, let's 
talk about why you'd want to calcu- 
late a square root. As an embedded 
programmer, I frequently use a square 
root to determine the root mean 
square (RMS) value of an AC voltage. 
This involves taking a square root. 

Other applications that use square 
roots are graphic algorithms and robot 
navigation systems. For example, to 
compute the distance between two 
points on a plane, you must use the 
Pythagorean theorem: 

'" d - yl&x* + Ay 2 ' 

So, why not just use the square root 
function in the G standard library? 
The sq rt ( ) C function computes the 
square root of a double-precision float- 
ing-point number, and returns a result 
of the same type. If all you need is the 
integer square root of a 16- or 32-bit 
integer, the overhead of the floating- 
point operations will slow down your 
application. For example, with the 
GCC compiler for AVR processors, 
sqrt( 50000) requires 2,254 cycles. 
The optimized 16-bit shift- subtract 
algorithm presented here takes only 
166 cycles, which is a 13x improvement. 

The wrong integer algorithm can also 
slow down your application. As demon- 
strated in this article, the difference 
between the slowest and fastest algo- 
rithms is an order of magnitude or more. 
Using commercial code doesn't guaran- 



tee an efficient algorithm. In fact, one 
of the commercial AVR compilers has 
integer square root functions that use 
an inefficient (though highly optimized) 
exhaustive-search algorithm. You might 
expect integer functions to be faster 
than floating point, but in my test of 
this compiler, taking the square root 
of 4,294,967,295 (the largest unsigned 
32-bit integer) took approximately 
655,000 cycles for the 32-bit integer 
function, 14,700 cycles for the floating- 
point function, and a mere 670 cycles 
for a 32-bit implementation of the 
shift-subtract algorithm presented 
later in this article. That's almost a 
l,000x difference between the slowest 
and fastest algorithms! 

The main topic in this article is 
integer square roots of integer argu- 
ments. All of the functions presented 
here take a 16-bit unsigned argument 
A and return a 16-bit unsigned integer 
result: x = VA 

Are you wondering why the return 
value is 16 bits instead of 8 bits? The 
width of the result is only half that of 
the argument, except for the routine 
that rounds the result instead of trun- 
cating. The square root of 65,281 (or 
higher) rounds to 256, which requires 
9 bits. In order to have a Common inter- 
face, I used only 16-bit return values. 

EXHAUSTIVE SEARCH 

I'll start with the exhaustive-search 
method, which is the simplest method. 
An exhaustive search basically tries 
every possible result until it finds the 
correct answer. Figure la tests succes- 
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Slues of x, starting at 1 and con- 
j to 255. For each loop, it com- 
x 1 with the argument A. If x 2 is 
liter than A, it terminates. (The test 
JL = handles the case, where x 
Ms over from 255 to because of the 



use of an 8-bit counter.) When the 
loop terminates, you know that the 
value of x is one greater than VA. The 
function returns x - 1 . 

Although simple, this algorithm has 
two drawbacks. First, it requires a 



uintl6_t' exhjnultiply 

(uintl6_t A) 
{ 

uint8_t x = 1; ■ 

while (((uintl6_t)x * x) 
<=- A) 



} 



++x; 
If (x 
{ 

break ; 

} 



- -X': . 
return (x); 



0) 



b) 

uintl6 t exh_sumodd 
(uint!6_t A) 



{ 



} 



uintl6_t x_sq = 1;' 
uintl6_t x_odd - 1; 

while (x_sq <= A) 
{ 

' x_odd += 2: 
x_sq += x_odd; 
if (x_sq = 0) 
{ 

break ; 



return (x odd » 1): 



o) 



uintl6_t exh_sumodd_opt 
(uintl6_t A) 

uintl6_t x_odd - 1; 

while (x_odd <= A) 
{ 

A — X_odd; 
x_odd +- 2; 

} . 

return (x odd » 1) ; 



multiplication every loop, which may or 
may not be a disadvantage, depending 
on whether or not your processor has a 
hardware multiply instruction. Second, 
like the other exhaustive-search algo- 
rithms, it loops times. For large val- 
ues of A, this can be significant, 
especially when the algorithm is 
extended to 32-bit arguments. 

Another exhaustive-search algo- 
rithm that eliminates the multipli- 
cation is the sum-of-odds method. 
In an exhaustive search, you need 
to know the value of x 2 for x = 1, 
2, 3, ... Because x is increasing by 1 
each loop, there's a shortcut to 
determine x 1 without performing a 
multiplication.; 

Suppose you know the value of x 1 
and you want to find the value of 
[x + If. Multiplying this out, you j 

(x + l) 2 =x J + 2x + l 



Figure 1— These functions use an exhaustive-search algorithm to calculate the square root of a 16-bit number. 
B&iaustive-search algorithms try all possible solutions until they find the correct solution. The simplest algorithm uses 
multiplication to compute the square of each trial (a). The sum-of-odds method 0mputes the square of each successive 
trial (b). An optimized sum-of-odds method improves performance '(c). Note that u1nt8_t and uintl6_t are type 
definitions for unsigned 8- and 16-bit integers. 



But 2x + 1 is just the index of suc- 
cessive odd numbers (i.e., first x + 1 
odd number). (2x gives successive 
even numbers. Adding 1 gives odd 
numbers.) 
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Figure lb shows the 
sum-of-odds algorithm in 
action. This function does- 
n't store the value of x. 
Instead it keeps track of x 2 
and the X th odd number (in 
x_sq and x_odd). After , 
each loop, it advances to 
the; next odd number by 
adding two. It then deter- 
mines x 2 by adding the odd 
number to the previous x 1 . 
Just like Figure la, it com- 
pares x 1 with A, but it does 
so without multiplication. 

After the loop termi- 
nates, the function needs 
to return x. But because it 
hasn't stored x in a vari- 
able, -it must calculate it. 
Because x_odd = 2x + 1, in 
order to convert x_odd to x, 
divide by 2 and truncate. 

Figure lc shows an opti- 
mized version of the sum- 
of-odds algorithm. Instead 
of keeping track separately of x 2 and 
comparing with A, this function sub- 
tracts the odd numbers directly from 
A. Because of this, the variable A is . 
always the difference between the origi- 
nal argument and x 2 . When that differ- 
ence goes negative, the loop terminates. 

As a final implementation of the 
sum-of-odds algorithm, I wrote an 
optimized function in AVR assembly 
for the Atmel ATmega8 processor. Not 
surprisingly, this optimized assembly 
code is the smallest and fastest of the 
exhaustive search implementations. 

BISECTION 

The sum-of-odds method is clever. 
After optimization, it has a fast inner 
loop, but it doesn't scale well to large 
arguments. In this section, I present 
methods that have much more consis- 
tent timing and that scale well: bisec- 
tion algorithms. 

Bisection continually divides the 
search space in half, producing 1 bit of 
the result for each iteration. Therefore, 
only eight iterations are required to gen- 
erate the square root of a 16-bit number, 
compared with up to 255 iterations for 
an exhaustive search. For 32-bit argu- 
ments, bisection requires 16 iterations, 
compared with up to 65,535 iterations 



a) 




b) 




c) 


uintl6 t bisect multiply 


uintl6 t bisect^shiftsub 


mntl6_t bisect_snsub_opt 




fnintlfi t AY 




fnintlfi t A1 


(mntl6_t A) 


r 

{ 




i 
V 




{ 
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uinto t Dit -uxou, 
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uintl6_t x_sh = 0; 




UjnXo l x , 




i n nf ft + y = fl ■ 
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uintl6_t bit_sh = x4000: 
uintl6_t trial ; 




while (bit 1-4) 
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x += bit : 
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while (bit sh 1-0) 
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{ 




if (((uintlb t)x * x; 




■f 
I 


trial = x_sh + bit_sh; 




> A) . 




trial - x_sq 


if (trial .<= A) 




{ . 




+ (x « (bnum + 1)) 


{ 




x -- bit ; 




+ (.1 « Kc * Dnum; ; ; 


A -= trial ; 




} 




if (trial <= A) 


x sh = trial + bit sh; 




bit »= 1 ; 




f 


} 




} 




x_sq = trial; 


x sh »= 1; 






x += (1 « bnum); 


bit sh »= 2; 




return (x) ; 




} 


} 


} 






-- bnum; 








} 

return (x); 


#if SQRT ROUND 

if (A > x sh) ++x sh: 






} 




#endif 








return (x sh): 

} 



Figure 2— 77iese functions use bisection to calculate the square root of a 16-bit number. Bisection algorithms divide the search 
space in half until they find the solution. You can use multiplication to compute the square of each trial {a). The shift-subtract method 
computes the square root of each successive trial (b). You can also use an optimized shift-subtract method for an improvement in 
performance (c). 



for an exhaustive search. A drawback to 
bisection is that it's more complex than 
an exhaustive search. 

As with the exhaustive-search algo- 
rithms, I begin with a simple algo- 
rithm that uses one multiplication per 
loop. I then replace the multiplication 
with more basic operations. I conclude 
with optimized implementations in C 
and assembly. 

Figure 2a shows bisection using one 
multiplication per loop. It starts by 
setting bit 7 of x to 1 (bit=0x80) and 
testing x 2 against A. If x 2 > A, it clears 
bit 7. It then does the same for bits 6 
through 0, and finally terminates after 
bi t becomes 0. During my tests, the 
algorithm's performance varied from 
poor to excellent. The Variations „ 
depended on whether or not the proces- 
sor supported hardware multiphcation. 



x 2 =(x p + 2 1 P 

x 2 =x 2 + 2x p 2 i + 2 21 

x 2 =x 2 + 2 i+1 x p + 2 2i 

X 2 = X 2 +[x p « (i + 1)] + (1« 34 ,;• 



Figure 3— A little algebra shows how the new value 
of x 2 can be calculated from the previous values, x„ 

and x 2 p . No multiplication is required, only addition 



The multiphcation hi Figure 2a can 
be removed. Each time through the 
loop, the function sets a single bit in x, 
which amounts to adding a power of 2. 
In Figure 2a, no variable indicates 
which power of 2. For now, assume 
variable i holds the power of 2 (the bit 
number) for each loop. During the 
first loop, i = 7. In subsequent loops, 
it's decremented down to 0. Each loop 
begins with x += bi t, which is equiv- 
alent to x += 2 1 . '"*.' -y-X 

Suppose you save the value of x 2 
from loop to loop. This will be handy 
in the next step, which calculates the 
new x 2 and compares it with A. Now 
get ready for a little bit of algebra. If x P 
is the previous value of x (so x p 2 is the 
previous value of x 2 ), you can calcu- 
late the new x 2 with the equations in 
Figure 3. ',4' 

That wasn't so bad, was it? After a 
few short steps, you've replaced the 
multiphcation in Figure 2a with left 
shifts and additions. All you need to 
do is save the previous value of x 2 and 
keep track of which bit is being set. 
Figure 2b is a direct implementation 
of this principle. It uses the variable 
bnum to keep track of the bit number 
(corresponding to i) and x_sq to save 
the value of x 2 . You can see that the • 
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computation of the intermediate 
trial variable exactly matches the 
final equation for x 2 . If the trial value 
doesn't exceed A, then the algorithm 
sets the bit in x and updates x 1 . 

Figure 2b is inefficient because of 
all the unnecessary shifting. Figure 2c 
is an optimized algorithm called the 
shift-subtract algorithm because the 
main operations shift and subtract. (For 
now, ignore the three lines starting with 
#i f SQRT_ROUND.),The first optimization 
is to eliminate the x_sq variable and 
subtract 2 i+1 x p + ^"from A. The result is 
that A always contains the error between 
the original argument and x 1 . If A is the 
original value passed to the function, 
then A - A - x 2 after each loop. 

The other optimization is to keep 
the preshifted values of 2 1+1 x p and V 1 
in the x_sri and bi t_sh variables (x 
shifted and bit shifted). Figure 4 
shows how tbis works. Even though 
there is no variable i in the function, 
each step is identified with an i value 
to make it easier to keep track. The 
bi t_sh variable represents 2 21 , so it's 
initialized to 0x4000 - 2 14 . The x_s h 
variable represents x shifted left by 
2 M . It's initialized to 0. 

Figure 4 shows only the bits of 
interest at each step in the calculation. 
A starts out with all bits unknown 
(represented by lowercase "a"), but as 
the algorithm progresses, some of the 
bits are forced to. 0. Take rny word for it. 
If you're skeptical, check out Figure 5 
for proof. 

The bi t_sh variable has a single 
1 bit, starting in bit 14, and shifted right 
by two positions each iteration. Finally, 
x vs h has 8 bits of interest that start in 
the high-order byte and are shifted right 
by one position for each iteration, ulti- 
mately ending in the low byte of the 
result. Each iteration generates 1 
bit of x_s h (represented by a low- 
ercase "x"). 

If you follow along in Figure 2c 
and Figure 4 simultaneously, you'll 
see that (x_sh + bit^sh) is ten- 
tatively subtracted from A in the 
comparison (tri a 1 <= A). . 
Remember that A has already had 
the previous value of x 2 subtracted 
from l^so this companson is effec- 
tively testing the original A g against 
the new tentative value of x 2 . If the 



i =7 
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i = 4 



i = 3 



1 = 2 
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Rgure 4—/ listed the variables [in the shift-subtract 
algorithm at thebeginning of each loop. Each iteration 
generates 1 bit of x and clears 1 bit of A. The i value is 
the bit number ofx that's being generated. Variable A 
contains the error after each iteration. 



comparison passes, then the algorithm 
subtracts the trial value from A and 
sets the bit in x^sh. 

VARIATIONS 

That completes my explanation of 
the shift-subtract algorithm, except for 
one variation. The algorithm as present- 
ed returns the truncated integer part of 
the square root. But what if you want to 
round the result instead of truncating . 
it? The basic algorithm would return 
49, which is V2,499. But 50 is a lot clos- 
er to the actual root. It turns but that 
you can use a simple test to round the 
result to the nearest integer. 
, , You must round up to the next inte- 
ger if the fractional part of the root is 
greater than 0.5. Round down other- 
wise; (It's impossible for the fraction 



A = aiii 


= A - (x- £i ) 2 = A -x 2 


+ 2xe, 




A - err. 


- A — A + Ie^Aq — e 2 = 


le^Ag 




A - en t 


<2(2'* 1 )2 8 - 2 2 *'* 1 ' 






A = err. 


< 2-io+ t 2 2 * i ^ 







Figure 5— The magnitude of A decreases as shown in Figure 3. err, is 
the value of the error (A) at each step as i decreases from 7 to 0. \ 

is the initial value of A, which must be less than 2" because it's a 
16-bit variable, x is the final result, so v? = A„, x -^A^ andi. < 2 s . 
Finally, e, is the error in x at Sach step, which is limited by me bits that 
have not been determined yet. Thus, e,<2 ; ". 



to be exactly 0-5 for the square root of 
an integer.) At the end of the shift- 
subtract algorithm, A = A^ - x 2 , which 
is greater than or equal to 0. If frac- 
tional part of VA b is greater than 0.5, 
then the following applies: 

A - (x + 0.5) 1 > 
A - x 2 - x - 0.25 > 
A > x + 0.25 

Because A and x are integers, the 0.25 
term can be ignored. The rounding 
test is simply to round up if A > x. In 
Figure 2c, you can see this implement- 
ed at the end of the function. To 
round the result, define SQRT_ROUND 
to be 1. To truncate it, make it 0. 

If you need to deal with fixed-point 
numbers, you might be wondering how 
the shift-subtract algorithm works with 
them. Fixed-point numbers have an 
implied binary point, so they can repre- 
sent numbers with an integer and a frac- 
tional part. They differ from floating- 
point numbers in that the position of the 
binary point doesn't vary at run-time. 

One way to handle fixed-point num- 
bers is to observe that the square root 
calculation divides the scale factor by 
two. So, if you prescale the input value 
to twice the desired output scale factor, 
the result will be correctly scaled. (This 
prescaling might require you to use a 
32-bit square root function for 16-bit 
numbers.) There are shift-subtract vari- 
ants that are tailored for fixed-point 
numbers, but I don't have space to 
cover those variants in this article. 

Another well-known algorithm for 
computing square roots is the 
Newton-Raphson algorithm (also 
known as the divide-and-average 
method, or Newton's method), which 
is the method I learned on the calcula- 
tor Way back when. I've seen pub- 
lished code that uses this method 
for integer square roots, but I rec- 
ommend against using it for this 
purpose for two reasons. 

First, it's relatively slow on 
embedded processors that lack a 
hardware divide instruction. If 
you look back at the shift-sub-. ; 
tract algorithm in Figure 2c, you 
can. see that it resembles the clas- 
sic binary division algorithm (in 
fact, the execution time is similar). 
The Newton-Raphson algorithm 
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Algorithm 1 


Instruction 
words ■ 


Cycles 


Normalized 


Average 


Maximum 


Average 


Maximum 


Exhaustive 
search 


SW Multiply 


54* 


20246 


31719: 


121.96 


.176.22 


HW Multiply 


16 


1716 


2558 


10.34 


14.21 


Sum-of-odds unoptimized 


19 


1688 


2819 


11.37 


15.66 


Sum-of-odds optimized C 


14 


1379 


2058 


8.31 


11.43 


Sum-of-odds optimized assembly 


n 


1208 


1802 


.7.28 


10.01 


Bisection 


SW Multiply 


42* 


1096 


1138 


6.60 


6.32 


HW Multiply 


16 


108. 


108 


0.65 


0.60 


Shift-subtract unoptimized 


63 


869 


952 


5.23 


5.29 


Shift-subtract optimized C 


26 


166 


180 


1.00 


1.00 


Shift-subtract optimized assembly 


25 


152 


166 


0.92 


0.92 



Table 1— Study the timing results for all of the 16-bit integer square root algorithms The results are normalized to 
the optimized. C language shift-subtract algorithm, f* SW multiply instruction words don't include the size of the 
library routine used ' to perform the multiplication (12 words for the AVR GCC library),) : 



requires several divide operations, 
each of which takes nearly the same 
amount of time as the entire shift-sub- 
tract calculation. 
; The other problem with the 
Newton-Raphson algorithm is that it 
doesn't always converge. It sometimes 
gives a result that's off by one. I rec^ 
ommend against the method for inte- 
ger square roots, although it can be a 
good algorithm for floating-point 
square roots. • 



RESULTS 

I compiled each of. the functions in 
Figures 1 and 2 using AYR GCC for the 
AlmegaS microcontroller, which sup- 
ports hardware multiply. I also compiled 
for the AT90S8515 microcontroller, 
which doesn't support hardware multi- 
ply, forcing the compiler to use software 
multiplication. Finally, I created hand- 
optimized AVR assembly code for the 
sum-of-odds and shift-subtract methods. 

Using a test program, I ran each func- 



tion with every argument from to 
65535. I also recorded the maximum and 
average number of clock cycles. (This 
took several minutes to run on a ■" 
14.7456-MHz Atmel ATmega8 micro- 
controller. ) The results are shown in 
Table 1. In addition to the number of 
cycles, I also included the number of 1-6- 
bit instruction words and the normalized 
number of cycles for each function. 
Normalization was achieved by dividing 
the cycles for each function by the cycles 
for the optimized C language implemen- 
tation of the shift-subtract algorithm. 
Different compilers or processors Would 
have had different results, but this test 
gives a good idea of relative performance. 

The first thing to note is that the 
worst bisection algorithm is faster than 
the best exhaustive-search algorithm. 
A fast inner loop executed 255 times 
doesn't beat a slower inner loop exe- 
cuted only eight times. 

The most surprising result is that the 
bisection algorithm with hardware multi- 
ply is the fastest of all. I really expected 
the optimized shift-subtract to be the 
fastest. But maybe this shouldn't be so 




only $98 my 1 




• 186 Processor @ 40 MH 

• DOS w/ Flash File System 

• Hardware Clock / Calendar 
•.512K.:Dl^^p^^pi 
r Console / Debug Serial Port 

^^P^H^^he pico - I /O 

|f|P2.P i'g ifalfi^rad^^^^^^ 

, C/C++* Quick Basic Drivers . 

-297-6073 Email sales@jkmicro.com 

at vy^^P^^M^c^^^^^S 




■■• (2) 16-bit Timers 

^^^^^^^^^^^^ 



J K microsystems 



Connecting the World of Electronic Design 



2006 



conference February 6-9, 2006 
exhibition February 7-8, 2006 

Santa Clara Convention Center 

Santa Clara, California 

www.designcon.com 



OFFICIAL SFOKSOR 




DIAMOND SP0NH0P 



International Engineering 
Consortium 

iuwiv.iec.org § 



64 Issue187 February 2006 



CIRCUIT CELLAR* 



www.clrcuitcellar.com 



surprising because the ATmega8 multi- 
plies in two cycles. (The shift-subtract 
algorithm takes six cycles just to do 
the shif,ting.) A processor with: a slow- 
er multiply instruction might give dif- 
: . ferent results. *>.f\j- V ^ ; «'* ; '- : v.- 

The 16-bit square root functions (in C) 
can be converted to 32 bits by doubling 
the widths of the argument, result, and 
intermediate values (i.e., change utnt8_t 
and uintl6_.t to ui ntl6_t and • 
uint32_t). I created 32-bit sum-of-odds 
and shift-subtract functions and measured 
the execution time for an argument of 
4,294,967,295 (the largest ui nt32_t). - 
The sum-of-odds algorithm took 917,518 
cycles. The shift-subtract algorithm 
took 670 cycles, a difference of more 
than three orders of magnitude! . _■ : 

So, now to the big question: Which 
algorithm is the best? 

Based on the timings in Table 1, the 
bisection algorithms are the best choice 
for general use. If I had to choose one inte- 
ger square root algorithm for my embed- 
ded tool kit, it would be the optimized 
shift-subtract algorithm in C because it , 
works well even on processors without 



hardware multiply. On a processor with 
hardware multiply, I'd consider the bisec- 
tion algonthm using a multiplication per 
loop, but only after testing the execution 
time versus the shift-subtract method. 
For the fastest possible algorithm for a 
particular processor, I'd compare opti- 
mized assembly language versions of the 
multiply and shift-subtract algorithms. 

The exhaustive-search algorithms did 
not perform well. Are there ever situa- 
tions in which an exhaustive search is 
appropriate? Yes, but such cases are 
uncommon. The optimized sum-of-odds 
algorithm is the most compact function. 
If code size is critical and execution time 
doesn't matter, the sum-of-odds algo- 
rithm might be acceptable. Or consider 
an application that computes the square 
rdot of slowly changing input data. In 
this sort of case, start each new search 
at the previous root, and do only a small 
' number of iterations per sample. You 
might find that a simple linear search 
outperforms bisection. 

I hope you find this article useful the 
next time you need an algorithm to com- 
pute an integer square root. Even if you 



don't need a square root algorithm, there 
are some useful conclusions to remember. 
The test results demonstrate that opti- 
mizing an inferior algorithm doesn't yield 
a superior solution. Sometimes what you 
think is an optimization is not. S 

Mike Dvoisky (bentbiker59-ccink@ 
yahoo.com) has B.S. and M.S. degrees 
in electrical engineering from the 
University of Missouri-Rplla. After 
working for seven years in VLSI 
design and verification, Mike moved 
to the dark side and became an 
embedded software engineer. He's 
currently an embedded software lead 
engineer at a Fortune 500 company. 
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To. download the code, go to ftp://ftp. 
circmtceUar.com/pub/Circuit_Gellar/ 
2006/187. 
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